function energy

%  comparison of energy using verlet and RK4
%  y'' = f(t,y)  with   y(0) = y0

% clear all previous variables and plots
clear *
clf

% parameters for calculation
n=8000;
t0=0; 
y0=[1 0];
tmax=4000;

%  calculate numerical solutions
t=linspace(t0,tmax,n);
h=t(2)-t(1);

%  RK4
[y_rk4,v_rk4]=rk4('de_f',t,y0,h,n);
for i=1:n
	e_rk4(i)=0.5*y_rk4(i)^2+0.5*v_rk4(i)^2;
end;

%  Verlet
[y_v,v_v]=verlet('de_f',t,y0,h,n);
for i=1:n
	e_v(i)=0.5*y_v(i)^2+0.5*v_v(i)^2;
end;

% get(gcf)
set(gcf,'Position', [1041 771 548 229]);

% plot energy
plot(t,e_v,'-r')
hold on
plot(t,e_rk4,'-b','LineWidth',1)
%plot(t,e_rk2,'-.k')
legend(' Verlet',' RK4','Location','East');

% define axes used in plot
axis([0 tmax 0 0.52])
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Energy','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
set(gca,'xtick',[0 1000 2000 3000 4000]);
box on
% Set the fontsize to 12 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off



% right-hand side of DE
function z=de_f(t,y)
z=[y(2) -y(1)];


% RK4 method
function [ypoints, vpoints] =rk4(f,t,y0,h,n)
y=y0;
ypoints=y0(1);
vpoints=y0(2);
for i=2:n
	k1=h*feval(f,t(i-1),y);
	k2=h*feval(f,t(i-1)+0.5*h,y+0.5*k1);
	k3=h*feval(f,t(i-1)+0.5*h,y+0.5*k2);
	k4=h*feval(f,t(i),y+k3);
	yy=y+(k1+2*k2+2*k3+k4)/6;
	ypoints=[ypoints, yy(1)];
	vpoints=[vpoints, yy(2)];
	y=yy;
end;

% RK2 method
function [ypoints, vpoints] =rk2(f,t,y0,h,n)
y=y0;
ypoints=y0(1);
vpoints=y0(2);
for i=2:n
	k1=h*feval(f,t(i-1),y);
	k2=h*feval(f,t(i),y+k1);
	yy=y+0.5*(k1+k2);
	ypoints=[ypoints, yy(1)];
	vpoints=[vpoints, yy(2)];
	y=yy;
end;


% Verlet method
function [ypoints, vpoints] =verlet(f,t,y0,h,n)
ypoints=y0(1);
vpoints=y0(2);
ya=y0(1);
va=y0(2);
for i=2:n
	y=ya+h*va-0.5*h*h*ya;
	v=va+0.5*h*(-y-ya);
	ya=y;  va=v;
	ypoints=[ypoints, y];
	vpoints=[vpoints, v];
end;